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A NEW EXTRAPOLATION CASCADIC MULTIGRID METHOD FOR 3D 
ELLIPTIC BOUNDARY VALUE PROBLEMS ON RECTANGULAR DOMAINS 


KEJIA PAN 1 * * * , DONGDONG HE § , AND HONGLING HU' 1 

Abstract. In this paper, we develop a new extrapolation cascadic multigrid method (ECMG/ C g), which makes 
it possible to solve 3D elliptic boundary value problems on rectangular domains of over 100 million unknowns 
on a desktop computer in minutes. First, by combining Richardson extrapolation and tri-quadratic Serendipity 
interpolation techniques, we introduce a new extrapolation formula to provide a good initial guess for the iterative 
solution on the next finer grid, which is a third order approximation to the finite element (FE) solution. And the 
resulting large sparse linear system from the FE discretization is then solved by the Jacobi-preconditioned Conjugate 
Gradient (JCG) method. Additionally, instead of performing a fixed number of iterations as used in the most of 
cascadic multigrid method (CMG) literature, a relative residual stopping criterion is used in our iterative solvers, 
which enables us to obtain conveniently the numerical solution with the desired accuracy. Moreover, a simple 
Richardson extrapolation is used to cheaply get a fourth order accurate solution on the entire fine grid from two 
second order accurate solutions on two different scale grids. Test results from three different problems with smooth 
and singular solutions are reported to show that ECMG ;C g has much better efficiency compared to the classical V- 
cycle and W-cycle multigrid methods. Since the initial guess for the iterative solution is a quite good approximation 
to the FE solution, numerical results show that only few number of iterations are required on the finest grid for 
ECMGy C g with an appropriate tolerance of the relative residual to achieve full second order accuracy, which is 
particularly important when solving large systems of equations and can greatly reduce the computational cost. It 
should be pointed out that when the tolerance becomes more cruel, ECMG 7C g still needs only few iterations to obtain 
fourth order extrapolated solution on each grid, except on the finest grid. Finally, we present the reason why our 
ECMG algorithms are so highly efficient for solving these elliptic problems. 
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1. Introduction. Elliptic boundary value problems arise in many areas of geophysical 
fluid dynamics. Consider the following model problem: 
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where a and /j are piecewise smooth functions on £1 and 0 < </3< /3 max for every x e £1, 

n is the outward unit normal to d£l, f : £1 —» go ■ Tzj —> 'X and gR : Fr —> % are assigned 
functions. Here £1 is some bounded rectangular domain in R 3 with Dirichlet boundary I /> and 
Robin boundary Yr. It is well known that the Neumann boundary condition corresponds to 
the extreme case, namely a - 0. 


*This research was supported by the Natural Science Foundation of China (Nos. 41474103, 41204082, 11402174 
and 11301176), the National High Technology Research and Development Program of China (No. 2014AA06A602) 
and the Natural Science Foundation of Hunan Province of China (No. 2015JJ3148). 

* School of Mathematics and Statistics, Central South University, Changsha 410083, PR. China (E-mail: panke- 
j ia @ hotmail. com) 

^Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223, 
USA 

§ School of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092, PR. China 
(E-mail: dongdonghe @ tongji.edu.cn) 

^Corresponding author (hhling625@163.com). College of Mathematics and Computer Science, Key Laboratory 
of High Performance Computing and Stochastic Information Processing (Ministry of Education of China), Hunan 
Normal University, Changsha 410081, PR. China 


1 






2 


K. J. PAN AND D. D. HE AND H. L. HU 


The elliptic boundary value problem can be approximated using different numerical tech¬ 
niques, such as finite difference (FD) and finite element (FE) methods. The resulting linear 
system can be solved efficiently using direct solver for problems with less than one million un¬ 
knowns. However, for 3D problems, even only a few hundreds grid points in each coordinate 
direction already leads to a system with millions of unknowns, one has to resort to an iterative 
method. For example, when solving direct current (DC) resistivity and electromagnetic mod¬ 
elling problems arising in geophysical applications, the grids designed to approximate huge 
realistic 3D geologies are usually enormous in order to represent correctly complex structures. 
Consequently, it is normally necessary to solve hundreds of millions unknowns in the forward 
problem. In addition, the forward problem has to be solved many times in the inversion of 
geophysical data [1, 2, 3], Therefore, it is critical that the 3D problem is solved very effi¬ 
ciently. The Multigrid technique [4, 5] is one of the most efficient strategies to solve the large 
linear system from discretized elliptic differential equations. The classical MG methods have 
been successfully applied to solve the Poisson equations [6, 7, 8, 9, 10, 11, 12], Helmholtz 
equation [13, 14, 15, 16, 17] and convection-diffusion equations [20, 22, 21, 18, 19]. How¬ 
ever, traditional MG methods (both geometric and algebraic) have to cycle between coarse 
and fine grids in order to accelerate the rate of convergence. Therefore, MG methods are 
difficult to implement in programming language. 

The Cascadic multigrid (CMG) method proposed by Deuflhard and Bornemann in [23] 
is a simpler multilevel method without coarse-grid correction. The CMG method uses CG 
solvers as the basic iteration methods on successively refined grids where the initial guesses 
are the linear interpolations of the approximate solutions on the previous grids. Nevertheless, 
the CMG method has the same optimal property compared to MG methods. Namely, the 
algorithm converges with a rate that is independent of the grid sizes and the numbers of grids 
levels [23, 24], In 2008, an extrapolation cascadic multigrid (ECMG) method was presented 
by Chen et al. in [42] for solving the second-order elliptic boundary value problems. This 
method proposes a new extrapolation formula to provide a better initial guess for the iterative 
solution on the next finer grid, which improves the convergence rate of the original CMG 
algorithm. However, as far as we know, the ECMG algorithm has mainly been used for 
solving the 2D elliptic boundary value problems in existing literature. But it is of more 
importance to solve the 3D problems efficiently and accurately arising in many engineering 
areas, such as geophysical exploration [ 1 , 2]. And it is nontrivial to extend the ECMG method 
from 2D to 3D. 

In this paper, we develop a new extrapolation cascadic multigrid method (ECMG JC g) for 
solving the 3D elliptic boundary value problems on rectangular domains. In our approach, 
the computational domain is discretized by a regular hexahedral grid, and a linear FE method 
is used to discretize the 3D elliptic problem. By extrapolating the numerical solutions on 
a coarse grid and a fine grid (with half the mesh size) to obtain a good initial guess of the 
iterative solution on the 8 vertices and 12 edge-midpoints of each coarse hexahedral element, 
which consists of 64 connected small hexahedral elements of the next finer grid (two times 
refined grid with one-fourth the mesh size), and then by using tri-quadratic Serendipity in¬ 
terpolation for the above 20 nodes to get the initial guesses on the other 105 (5 3 - 20) nodes 
of such coarse hexahedral element, we are able to obtain a third order approximation to the 
FE solution on the next finer grid. And the resulting large sparse linear system from the FE 
discretization is then solved by the Jacobi-preconditioned Conjugate Gradient (JCG) solver 
using the obtained initial guess. Additionally, a tolerance related to relative residual is in¬ 
troduced in the JCG solver. Moreover, a simple Richardson extrapolation is used to obtain 
cheaply a fourth order accurate solution on the entire fine grid from two second order accurate 
solutions on two different scale grids (current fine and previous coarser grids). Finally, our 
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method has been used to solve 3D elliptic problems with more than 135 million unknowns 
in about 1 minute on a small server with 32GB RAM installed. Compared to the existing 
ECMG method, our new ECMG approach is advantageous due to the following reasons: 

• Instead of performing a fixed number of iterations as used in the usual CMG meth¬ 
ods, a relative residual tolerance is used in our iterative solvers, which enables us 
to not only avoid the difficulty of choosing the number of iterations on each grid, 
but also allows us to obtain conveniently the numerical solution with the desired 
accuracy. 

• The second important point is to employ JCG method as MG smoother, which re¬ 
quires only a little additional work than CG, while JCG smoother typically yields 
better convergence properties than CG smoother. 

• A fourth order extrapolated solution on the entire fine grid is constructed to greatly 
enhance the accuracy of the numerical solution. 

• we can clearly explain why only few number of iterations are required on the finest 
grid to achieve full second order accuracy for our method through defining a ratio. 

The rest of the paper is organized as follows: Section 2 gives the description of the FE 
discretization for the 3D elliptic boundary value problem. Section 3 reviews the classical 
V-cycle and W-cycle MG methods. In Section 4, we first present some extrapolation and 
interpolation formulas, and then develop a new ECMG method to solve 3D elliptic boundary 
value problems. Section 5 contains the numerical results to demonstrate the high efficiency 
and accuracy of the proposed method. And conclusions are given in the final Section. 

2. Finite element discretizations. In this Section we introduce a linear FE discretiza¬ 
tion to (1.1) which we shall use in the construction of MG methods. For simplicity, we assume 
that eq.(l.l) has homogeneous Dirichlet boundary conditions, namely go = 0. We introduce 
a bilinear form, a(-, •): x H l D (Q) —> R , defined in the usual way 

(2.1) a(u,v)= I /3Wu-Wvdx + I auvdx, 

J Q 

where = {v e //'(£!) : v| ro = 0}, and v|r D = 0 are in the sense of trace. 

Then the weak formulation of the problem (1.1) is to find u e /I],(£!) such that 

(2.2) a(u, v ) = /(v), Vv e H l D (Cl), 
where 

(2.3) f(v) = f fvdx + f g R vds. 

Assuming that O is partitioned by an uniform hexahedral mesh T with characteristic 
mesh size h, namely Q = {J reTli Then a piecewise linear FE space, 14 c H(Cl) can be 
defined by 

(2.4) 14 = {v e ff*,(£2): v| T 6 Pi(t), Vr e T h ), 

where P\ denotes the set of linear polynomials. 

Let {0,}^ be the standard nodal basis functions of FE space 14 and writing the FE 
solution as the linear combination of the basis functions as: u h = Yi'Li u h,i<t>i , we can obtain 
the following algebraic systems for the variational problem (2.2): 


( 2 . 5 ) 


Ah^h — fhi 
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where 

(2.6) A h = (a(0i, <pjj ), u h = (u hJ ), f h = if {(pi)). 

In this paper, we focus on how to solve the FE equation (2.5) efficiently using a new extrapo¬ 
lation cascadic multigrid method. 

3. Classical multigrid method. The MG method is very useful in increasing the effi¬ 
ciency of iterative methods (such as Gauss-Seidel or weighted lacobi) for solving the large 
systems of algebraic equations resulting from FE discretizations of elliptic boundary value 
problems. The main idea behind the MG methods is to split the errors into high frequency 
errors and low frequency errors and solve them in different subspaces. To be more specific, 
the high frequency errors are solved on the fine grid, while the low frequency components are 
solved on the coarse grid. A good introductory text on multigrid is the book by Briggs [4]; 
and a more advanced treatment is given by Trottenberg in [5]. 


Algorithm 1 Classical MG Method (Recursive Definition): Uh <= MGMM/,, Uh,fh) 
1: if h == //(Coarsest grid) then 
2: uh <= DSOLVE(A//« w = f H ) 

3: else 

4: Perform N pre pre-smoothing iterations with the current estimate m/, 

f/i — fh A/jUh 

6 : r 2h = Rl h r h 

7 : e 2 h = 0 

8: for i = 1 to N cycles do 

9: e 2 h MGM(A 2 /,, e 2 h , r 2h ) 

10 : end for 

11: Cji R oj^lh 

12 : Ufj = U h + e h 

13: Perform N pos , post-smoothing iterations with initial guess «/, 

14: end if 


Algorithm 2 Multigrid solver for the finest grid system A/,iq, = //, 

1: M/, = 0 

2: while || A h u h - f h \\ 2 > e\\f h \\ 2 do 
3 : Uh <= MGM(A/j, m / j , y/ 2 ) 

4 : end while 


The MG method can be defined recursively as follows. On the coarsest grid, a direct 
solver DSOLVE is used since the size of the linear system is small, see line 2 in the Algorithm 
1. The procedure MGMtA/,,«/,,/;,) can be applied to the fine grid system (2.5) in the solution 
of which we are interested. The MG iteration with procedure MGM is repeated on the fine 
grid until the fine residual is considered small enough, see Algorithm 2 for details. 

To fully specify the Algorithm 1, one needs to specify the numbers of pre- and post¬ 
smoothing steps on each grid N pre and N post , the restriction (fine-to-coarse) operator/?^', the 
prolongation (coarse-to-fine) operation P h lh , and the recursion parameter N cyc i es . One iteration 
of a MG method, from the finest grid to the coarsest grid and back to the finest grid again, 
is called a cycle. The exact structure of MG cycle depends on the parameter N cyc i es , which 
is the number of two-grid iterations at each intermediate stage, greater or equal to one. With 
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Grid level 



2 

: 


V-cycle W-cycle CMG EXCMG 

Fig. 3.1: The four level structure of the V- and W-cycles, CMG and ECMG methods. In the 
diagram, • denotes pre-smoothing steps, o denotes post-smoothing steps, f denotes prolon¬ 
gation (usually defined by linear interpolation), X denotes restriction, ff denotes extrapolation 
and quadratic interpolation, and ■ denotes direct solver. 


Ncycles = 1, the so-called V-cycle is generated, while N cyc i es - 2 leads to W-cycle. The four 
level structures of the V- and W-cycles are illustrated in Fig. 3.1 

The proper choices of both residual restriction and error prolongation operators play a 
key role in the successful development of the MG method. For prolongation operator P^ h , we 
use the standard tri-linear interpolation operators which preserves the symmetry of the true 
solution exactly. And the 27-point fully weighted residual restriction is used in this study, 
which is a natural extension from the 2D 9-point full weighting restriction to 3D space. The 
restriction operator Rj t h is defined as follows [ 5 ]: 


+i 

(3.1) Rfu h {x,y,z) = Z WijkUhix + /Ax, y + jAy,z + kAz), 

i,j,k =-1 

where 


1/8 

if 

i = j = k = 0, 

1/16 

if 

\i\ + I/I + 1*1 = 1, 

1/32 

if 

lil + I/I + 1*1 = 2, 

1/64 

if 

lil + I/I + 1*1 = 3. 


In this manner, the weighting satisfied the conservation property of the integrals. 

4. Extrapolation cascadic multigrid methods. The CMG method proposed by Deufl- 
hard and Bornemann in [ 23 ] is a type of method which requires no coarse grid correction at all 
that may be viewed as a one-way MG (see Figure 3.1). Because no coarse grid correction is 
needed, computational cost is greatly reduced compared to the classical MG methods. Since 
the 1990s, the method received quite a bit of attention from researchers because of its high 
efficiency and simplicity [ 24 , 25 , 26 , 27 , 28 , 29 , 30 , 31 , 32 , 33 , 34 , 35 , 36 , 37 , 38 , 39 , 40 , 41 ], 
By using Richardson extrapolation and bi-quadratic interpolation techniques, an extrapola¬ 
tion cascadic multigrid (ECMG) method for elliptic problems was first proposed in 2008 by 
Chen et al. in [ 42 , 43 ]. And some 2D numerical experiments are presented in [ 44 ] to show 
the advantage of ECMG in comparison with CMG method. 
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4.1. Description of the ECMG algorithm. The key ingredients of the ECMG method 
are extrapolation and quadratic interpolation (see Figure 3.1), which are used to provide a 
better initial guess of the iterative solution on the next finer grid than one obtained by using 
linear interpolation in CMG. The ECMG algorithm is given as follows [44]: 


Algorithm 3 ECMG : itj, <= ECMG(A/j, f ,, L, mj) 

1: u H <= DSOLVE(A//M// = f H ) 

2 : Uh/2 <= DSOLVE(A h /2 M///2 = fm 2 ) 

3 : h = H/2 
4: for i - 1 to L do 
5: h = h/2 

6: Uh — EXPy/ tt // e (i/2/j, U 4 h) 

7: for j = 1 to nij do 

8: U / 7 4= CG(A/ Z , M/j, fh) 

9: end for 

10: end for 


In Algorithm 3, procedure CG denotes conjugate gradient iteration, and procedure EXP 
denotes extrapolation and bi-quadratic interpolation operator on three levels of the embedded 
rectangular grids, which is used to get a good initial guess of the iterative solution on the next 
finer grid, see line 6 and line 8 in the Algorithm 3. The positive integer L, total number of 
grids except first two embedded grids, indicates that the mesh size of the finest grid is ^r. 
And on the first two coarse grids, a direct solver DSOLVE is used since the size of the linear 
system is small, see line 1-2 in the Algorithm 3. 

The number of iterations on each grid is defined by 

(4.1) nij — \iriL -/? L ; ], 7 = 2,3,■ • • ,L, 

where m.L is the iteration count on the finest grid, fl is usually set to be between 2 and 2 d , 
and [jc] denotes the smallest integer greater than or equal to x. In [53], authors compared 
different settings for mr and /i, and found that ECMG with nij - 4 x 4 L ~ J can achieve high 
accuracy with quick convergence. In 2007, Shi, Xu and Huang [37] proposed an economical 
CMG method, which uses a new criteria for choosing the smoothing steps on each level, when 
d - 2 , that is 

(i) if j > L () , then nij = \m L ■ /3 L ~ j ~\, 

(ii) if./' < To, then nij — \(L - (2 - e 0 )j)h^ 2 l 

Here 0 < eo < 1 is a fixed positive number, and //!/ = mo(L - To) 2 . Note that in standard 
cascadic multigird algorithm //// = nioL 2 . The level parameter Lq, which depends on L, mi,/3 
and the size of the coarest level H, can be defined as the largest integer satisfying 


.. .. [ Tlog/3 + log m L + 41og// L \ 

(4.2) To ^ mini-, — 

I log/? + 4 log 2 ’2j 

To our best knowledge, the ECMG method has mainly been applied to solve 2D elliptic 
boundary value problems, but it is very important to solve the 3D elliptic problems efficiently 
and accurately in many engineering areas and the extension of ECMG method from 2D to 
3D is nontrivial. The advantage of the ECMG method will be much more reflected when 
it is applied to 3D problems since the size of linear discrete system is much larger. Next, 
we propose a new ECMG method for solving 3D elliptic problems which is stated in the 
following algorithm 4: 
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Algorithm 4 New ECMG : ( Uh , Uh ) <= ECMG(A/,, //,, L, e) 
l: u H <= DSOLVE(A//M// = f H ) 

2: Uh /2 <= DSOLVE(A///2Mh/2 = ffi/l) 

3 : h = H /2 
4: for i = 1 to L do 
5 : h = h/2 

6: Wh = EXP/,„ i(< ,(M 2 /i, M 4 a) > w/i = w/, is used as the initial guess for JCG solver 

7: while \\A h u h - f h \\ 2 > e ■ \\f ,|| 2 do 

8 : M/ 7 4= JCG(A/ ? , Uh, fil) 

9: end while 

10: S/j = EXP i rll e{uh, « 2 /i) > S/, is a higher-order approximate solution 

11 : end for 


In Algorithm 4, EXP/,„„ ( ,(m 2 a, « 4 /,) denotes a third order approximation to the FE solu¬ 
tion ui, obtained by Richardson extrapolation and tri-quadratic Serendipity interpolation tech¬ 
niques from the numerical solutions on embedded hexahedral grids, while EXP„ Mf .(M/,, u 2 h) 
denotes a fourth order extrapolated solution on fine grid with mesh size h from two second 
order numerical solutions m/, and u 2 i,. The detail procedures of extrapolation and Serendip¬ 
ity interpolation on embedded hexahedral grids are described in the next Section 4.2. The 
differences between our new ECMG and the existing ECMG are illustrated as follows: 

(1) Instead of performing a fixed number of CG iterations (see line 7 in Algorithm 3) as 
used in the usual CMG methods, a tolerance e related to the relative residual is used 
in our iterative solvers (see line 7 in Algorithm 4), which enables us to not only avoid 
the difficulty of choosing the number of iterations on each grid im, but also allows us to 
conveniently obtain the numerical solution with the desired accuracy. 

(2) Algorithm 4 takes JCG as the MG smoother, which is written as the line 8 in the Al¬ 
gorithm 4. The Jacobi preconditioner, which consists only of the diagonal elements of 
the matrix, can be constructed from the coefficient matrix without any more extra work. 
Although JCG requires a little additional work than CG, JCG smoother typically yields 
better convergence properties than CG smoother (see Table 5.1- 5.3 in Section 5). 

(3) A fourth order extrapolated solution Uh (see Table 5.6- 5.9 for details) on the current 
grid based on numerical solutions of two level grids, current grid and previous (coarser) 
grid, is constructed in Algorithm 4 (see line 10) in order to enhance the accuracy of the 
numerical solution Uh- 

In the following, we denote the above new ECMG method as ECMG /a , and denote 
ECMG C g for the Algorithm 4 where the line 8 is replaced by 4= CG (A/,, u/,, fk). 

4.2. Extrapolation and quadratic interpolation. It is well known that the extrapo¬ 
lation method, which was established by Richardson in 1926, is an efficient procedure for 
increasing the solution accuracy of many problems in numerical analysis. In 1983, Marchuk 
and Shaidurov [ 45 ] studied systematically the application of this method in finite difference 
method. Since then, this technique has been well demonstrated in the frame of the FE method 
[ 46 , 47 , 48 , 49 , 50 , 51 ], 

In this Section, we will explain how to use extrapolation and interpolation techniques 
to obtain a high order accurate approximate solution to the problem (1.1), and a high order 
accurate approximation to the FE solution. And it can be regarded as another important 
application of the extrapolation method to provide a good initial guess of iterative solution. 
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4.2.1. Extrapolation for the true solution. For simplicity, let us first consider the three 
levels of embedded grids Z,-(i = 0,1,2) with mesh-sizes /;, = /z 0 /2' in one dimension. Let 
e' - U l - u be the error of the linear FE solution U' e V/,., it can be proved that under certain 
smoothness conditions the error at node x k has the the following asymptotic expansion 

(4.3) e‘(x k ) = A(x k )h 2 + 0(h% 

where A(x) is a suitably smooth function independent of /;,. 

It is well known that the extrapolation methods can only offer a fourth-order accurate ap¬ 
proximation to the true solution at coarse grid points. The Richardson extrapolation formula 
at coarse grid points { xj , xj +\) can be written as 

_ 4 jj l - jj° 

(4.4) U l k :- * 3 * = u(x k) + O(h 4 0 ), k = j,j + 1. 

In fact, by using the linear interpolation formula, one can also obtain a fourth-order ac¬ 
curate approximation at fine grid points. Chen and Lin [49] proposed a fine grid extrapolation 
formula 

(4.5) U) +l/2 := t /] +1/2 + l -(U) - U°j + U) +l - t/° +1 ) - u(x j+1/2 ) + O(h 4 0 ), 

which can give directly the higher order accuracy at fine grid points. 

From (4.3), we can obtain 

(4.6) A(x k ) = -l(t/° - Ul) + 0(h 2 Q ), k = j,j + 1. 

3n 0 

Applying the error estimate of the linear interpolation, 

(4.7) A(Xj+ 1 / 2 ) = —(A(xj) + A(xj+ 1 )) + O(Iiq). 

Substituting eq. (4.6) into eq. (4.7) yields 

(4.8) A(Xj + \/2) = ^AU) - U)) + -^(t /« +1 - U) +1 ) + O(h 2 0 ). 

3 n o ^ n o 


Since 

( 4 . 9 ) U j+i/2 = U ( X J+ 1 / 2 ) + ^ A (x j+ i/2)hl + 0 (/ iq ), 

by using (4.8), the extrapolation formula (4.5) is obtained. 

4.2.2. Extrapolation for the FE solution. Next, we will explain, given the FE solutions 
U° and U 1 , how to use the extrapolation method to obtain a third order (to be proved in 
subsection 4.4) approximation W 2 for the FE solution U 2 rather than the true solution u. 

Adding one midpoint and two four equal division points, the coarse mesh element (xj, x/+i) 
is uniformly refined into four elements as shown in Fig. 4.1. And we obtain a five points set 

jx 2 , Xy+i/4, Xj.y | /2, Xj+ 3/4, Xj.\. ] J 


belonging to fine mesh Z 2 . 
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u 1 


u,\ 


U 1 , 


u 


u 


j+1/4 


u: 


j+1/2 


u: 


j+3/4 


u: 


j+i 


j j+1/4 j+1/2 j+3/4 j+1 

Fig. 4.1: Three embedded grids in ID. 


Assume there exists a constant c such that 


(4.10) ell 1 + (1 - c)U° - U 2 + O(Iiq), 

that is to say, we use a linear combination of U° and U 1 to approximate the FE solution 
U 1 . Substituting the asymptotic error expansion (4.3) into (4.10), we obtain c = 5/4 and an 
extrapolation formula 

5 f/ 1 _ u° 

(4.11) W 2 := * 4 k = U\ + 0(h 4 ), k = j, j + 1, 

at nodes Xj and Xj+ \. To derive the extrapolation formula at mid-point x ;+ i/ 2 , we have 
(4-12) U 2 j+m = U) +1/2 - ^ A(x j+l/2 )hl + O(h 4 0 ), 


by using equation (4.3). Substituting eq.(4.8) into eq.(4.12) yields the following mid-point 
extrapolation formula. 


(4.13) W) +m := U) +m + -(U) - U° + U) +l - U° j+1 ) = U 2 j+l/2 + 0(h 4 ). 


Once the three initial values W 2 , Wj +1 and VF " +1 ^ 2 are obtained, we can get the following four 
equal division point extrapolation formulas by using quadratic interpolation, 


<"■ 14 > Km : = 4[(9i/J + 12(3 U° + (/;,,)], 

(4.15) Kw ■= ^K w ;+I + * u °>]' 


4.3. Three dimensional case. In this subsection, we explain how to obtain a third order 
accurate approximation W 2 to the FE solution U 1 , and a fourth order accurate approximate 
solution U l to the problem (1.1) for embedded hexahedral grids as shown in Fig. 4.2. 

The construction process of the approximation W 2 are as follows: 

Corner Nodes (1,5, 21,25,101,105,121,125): The approximate values at 8 corner nodes 
can be obtained by using the extrapolation formulae (4. 1 1) 
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Fig. 4.2: Three embedded grids in 3D. 


Midpoints in ^-direction (3, 23,103,123): The approximate values at these 4 midpoints 
can be obtained by using the mid-point extrapolation formulae (4.13) in x- 
direction. 

Midpoints in y-direction (11,15, 111, 115): The approximate values at these 4 midpoints 
can be obtained by using the mid-point extrapolation formulae (4.13) in y- 
direction. 

Midpoints in ^-direction (51, 55, 71, 75): The approximate values at these 4 midpoints 

can be obtained by using the mid-point extrapolation formulae (4.13) in ^-direction. 

Other 105 points: The approximate values of remaining 105 (5 3 - 20) grid points can be 
obtained by using tri-quadratic Serendipity interpolation with the known 20-node (8 
corner nodes and 12 mid-side nodes) values. 

The tri-quadratic Serendipity interpolation function in terms of natural coordinates (£, q, if) 


is 


(4.16) 



where the shape functions /V, can be written as follows [52]: 
(4.17) 



i = 1,5,21,25,101,105,121,125 


i(W 2 )(l+/7,77)(l+£a 


i = 3,23,103,123 


Ni(g,ri,0 = \ , 

-(l-ift(l+ 60 (l+&). 


i= 11,15,111,115 


i(l-f 2 )(l+6f)(l + 77,7,), 


i = 51,55,71,75 


where (f, rji, £,) is the natural coordinate of node i. For example. 


(4.18) 



Remark 4.1 . Extrapolation and tri-quadratic interpolation defined by Eq.(4. 11), Eq. (4.13) 
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and Eq.(4.16) are local operations, which can be be done very cheaply and effectively. 

Remark 4.2. Since the center of each face such as the point labeled 53 is the midpoint 
of two face diagonals, we can also use the mid-point extrapolation formulae (4.13) to obtain 
two different extrapolation values, and take their arithmetic mean as the approximated value. 
A similar procedure can be done for the center point of the cube labeled 63, which can 
be viewed as the midpoint of space diagonals. Then the approximated values of remaining 
98 (5 i - 27) points can be obtained by using tri-quadratic Lagrange interpolation with the 
known 27-node values. This tri-quadratic Lagrange interpolation also gives a third order 
approximation of the FE solution on the finer grid followed by a similar error analysis in the 
next subsection, in addition, numerical results show that it has almost the same accuracy as 
the tri-quadratic Serendipity interpolation with the known 20 nodes. Thus, we only present 
the numerical results by using the method of tri-quadratic Serendipity interpolation in the 
next Section. 

When constructing the fourth order accurate approximate solution f/ 1 based on two sec¬ 
ond order accurate solutions U° and U 1 , the Richardson extrapolation formula (4.4) can be 
directly used for 8 coarse grid points, while the fourth order extrapolation formula (4.5) can 
be used for the 12 centers of edges, the 6 centers of faces as well as the center point of the 
coarse hexahedral element. 

Remark 4.3. In the existing literature, there are already some extrapolation related meth¬ 
ods which use the two numerical solutions from the tw’o level of grids to obtain high order 
accuracy solution on the fine grid. For example, Wang et al. [11 ] first used the Richardson 
extrapolation for two fourth order compact finite difference solutions on two level of grids 
to achieve sixth order accuracy on the even number of grids of the fine grid, however, the 
method of [11] relies on an operator based interpolation iterative strategy for odd number 
of grids in order to achieve the sixth order accuracy on the whole fine grid, which increases 
the computational cost. Our extrapolation formulas (4.4)-(4.5) can be used to obtain fourth 
order accurate solution on the whole fine grid from two second order accurate solutions on 
two level of grids directly and cheaply without any extra work. 

4.4. The error analysis of initial guess W 2 . Let e - W 2 - U 2 be the approximation 
error of extrapolation formulas, and assume that it has a continuous derivative up to order 3 
on interval [xj,Xj+ 1 ]. From (4.11) and (4.13) we obtain the equation 


e(x k ) = O(h 4 0 ), k = j,j + 1 / 2,7 + 1. 


(4.19) 


From polynomial interpolation theory, the error of quadratic interpolation ff can be repre¬ 
sented as, 



(4.20) 


where f is a point of [x ; , Xj+ 1 ] that depends on x. Especially at four equal division points we 
have 



and 



It follows from eq. (4.21), eq.(4.22) and eq.(4.19) that 

(4.23) e(x k ) = I 2 e(x k ) + R 2 (x k ) = O(h 3 0 ), k = j + 1/4,;' + 3/4, 
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Table 5.1: Comparison of the numbers of iterations and CPU times (in seconds) of four 
different MG methods with e = 10” 8 for Problem 1. 


Mesh 

ECMGj Cg 

ECMG cg 

V(U) 1 

W(2,l) 2 

32 x 32 x 32 

7 

58 

26 

432 

64 x 64 x 64 

10 

82 

26 

216 

128 x 128 x 128 

18 

93 

26 

108 

256 x 256 x 256 

3 

58 

26 

54 

512x512x512 

3 

9 

26 

27 

CPU (s) 

147 

169 

452 

466 


1 The number of MG V(l,l) cycle is 13. 

2 The number of MG W(2,1) cycle is 9. 


Table 5.2: Comparison of the numbers of iterations and CPU times (in seconds) of four 
different MG methods with e = 10 9 for Problem 1. 


Mesh 

ECMGj Cg 

ECMG cg 

V(U) 1 

W(2,l) 2 

32 x 32 x 32 

8 

71 

30 

480 

64 x 64 x 64 

9 

98 

30 

240 

128 x 128 x 128 

16 

150 

30 

120 

256 x 256 x 256 

78 

162 

30 

60 

512x512x512 

3 

50 

30 

30 

CPU (s) 

162 

284 

489 

545 


1 The number of MG Vf 1,1) cycle is 15. 

2 The number of MG W(2,1) cycle is 10. 


which means that four equal division point extrapolation formulas (4.14) and (4.15) are only 
third-order approximations. 

The above error analysis can be directly extended to 3D case (also see numerical verifi¬ 
cation of Table 5.4-5.10 in the next Section). In addition, eq.(4.22) implies that the error e(x) 
forms a high-frequency oscillation in whole domain, however, it can be smoothed out after 
very few iterations (see Fig. 5.1 for details). 

5. Numerical experiments. In this Section, we will illustrate the efficiency of ECMG 7 - cg 
by comparing to ECMG cg as well as the classical V-cycle and W-cycle MG methods, and 
present numerical results for three examples with smooth and singular solutions obtained 
by the most efficient algorithm ECMG ;r? . In our experiments, we use the most popular 
smoothing method Gauss-Seidel smoother as the relaxation smoother of classical MG meth¬ 
ods, which usually leads to a good convergence rate. Our code is written in Fortran 90 and 
compiled with Intel Visual Fortran Compiler XE 12.1 compiler. All programs are carried out 
on a server with Intel(R) Xeon(R) CPU E5-2680 (2.80GHz) and 32G RAM. 

Problem 1. The test Problem 1 can be written as 


(5.1) 


d 2 u d 2 u d 2 u 3 7 . 7r 7r n 

+ 3f + a? = -r “*^>s,n(-,)sm(- z ). 


ill ii = [0, l] 3 . 


where the boundary conditions are 


|^(1 ,y,z) = |^(x, l,z) = 1) = 0. 

on on on 


(5.2) u(0,y,z) = u(x,0,z) = u(x,y, 0) = 0, 
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Table 5.3: Comparison of the numbers of iterations and CPU times (in seconds) of four 
different MG methods with e = 10 10 for Problem 1. 


Mesh 

ECMGy cg 

ECMG cg 

V(U) 1 

W(2,l)“ 

32 x 32 x 32 

9 

76 

32 

528 

64 x 64 x 64 

9 

124 

32 

264 

128 x 128 x 128 

26 

172 

32 

132 

256 x 256 x 256 

112 

258 

32 

66 

512x512x512 

3 

80 

32 

33 

CPU (s) 

177 

371 

522 

534 


1 The number of MG V(l,l) cycle is 16. 

2 The number of MG W(2,1) cycle is 11. 


The analytic solution of eq. (5.1) is 

7T 7T 7T 

u(x, y, z ) = stn( -x) stn( -y) sin( -z), 

where the exact solution is a sufficiently smooth function which has arbitrary order smooth 
derivatives and the variation of the function is the same in three directions. 

Using 7 embedded grids with the coarsest grid 8 x 8 x 8, we give the CPU time and the 
number of iterations at each level of grid for ECMG 7 - cg , ECMG cg , the classical V-cycle and 
W-cycle MG methods using three different tolerances e = 10 s , 10~ 9 and 10 1,1 in Table 5.1, 
Table 5.2 and Table 5.3, respectively. Here the CPU time is total computational time which 
not only includes the time to solve the linear system by MG methods but also includes all 
other necessary computational time in the program. Since the first two coarse level of grids 
are solved by direct solver, we only list the results start from the third level of grid 32x32x32. 

As we can see from Table 5.1 to Table 5.3 that the ECMG ;cg uses the minimum number 
of iterations and minimum CPU time and has the best efficiency among all MG methods, 
ECMG cg is a little worse than ECMG Jcg , while the classical MG with W-cycle has the worst 
efficiency and has the similar efficiency as the classical MG with V-cycle. In addition, we 
can see that when the control parameter e becomes smaller, the number of iterations and the 
CPU time for all MG methods increase as expected, and the efficiency of ECMGj Cg becomes 
much better than ECMG cg while the classical MG with V-cycle still has the similar efficiency 
as the classical MG with W-cycle. It should be noted that when the tolerance becomes 10 l0 , 
only 3 JCG iterations are required on the finest grid 512x512x512 (more than 135 million 
unknowns) for ECMG ycg method, which costs less than 3 minutes (see second column of Ta¬ 
ble 5.3). Moreover, numerical comparisons for efficiency are also carried out for the next two 
examples (not listed in the paper), results show that ECMG 7 - cg always has the best efficiency 
among these four MG methods. 

For the sake of clarity, however, it should be mentioned that the process of solving the 
linear system only required a small portion of the total computing time compared to the effort 
spent in discretizing the problem, i.e., computing the stiffness matrices and load vectors. 
Additionally, since we are using different MG methods to solve the same linear system, the 
accuracy of all methods should be similar under the same tolerance. Thus, in the following 
part of this paper, we present only the numerical results obtained by ECMGj Cg since it is the 
most efficient method. 

We present the numerical results for Problem 1 obtained by ECMG jr ,, with e = 10 8 
in Table 5.4-5.5, and with e - 10 9 in Table 5.6-5. 7, where “Iters” denotes the number of 
iterations needed for the JCG solver to achieve that the relative residual is less than the given 
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Table 5.4: Errors and convergence rates with e - 10 8 in Li norm for Problem 1. 


Mesh 

Iters 

RRe 

Ilf//, - 

"lb 

Ilf//, - 

m|| 2 

\\W h - 

UhWi 

Hi 

Error 

Order 

Error 

Order 

Error 

Order 

32 x 32 x 32 

7 

9.99(-9) 

1.42(-4) 


1.96(—7) 


2.54(-5) 


1.79(-1) 

64 x 64 x 64 

10 

8.14(—9) 

3.55(-5) 

2.00 

1.24(-8) 

3.98 

3.18(—6) 

2.99 

8.96(-2) 

128 x 128 x 128 

18 

7.68(-9) 

8.87(-6) 

2.00 

7.83(—10) 

3.99 

3.99(—7) 

3.00 

4.50(-2) 

256 x 256 x 256 

3 

6.26(-9) 

2.22(-6) 

2.00 

4.74(-10) 

0.72 

4.99(-8) 

3.00 

2.25(-2) 

512x512x512 

3 

6.00(-9) 

5.55(-7) 

2.00 

4.69(-10) 

0.02 

6.25(-9) 

3.00 

1.13(—2) 


Table 5.5: Errors and convergence rates with e = 10 8 in norm for Problem 1. 


Mesh 

Iters 

RRe 

Ilf//, - 

W| |cX3 

Ilf//, - 

|co 

Wh - 

UhWoc 

Error 

Order 

Error 

Order 

Error 

Order 

32 x 32 x 32 

7 

9.99(-9) 

4.02(-4) 


l.ll(-6) 


2.54(-5) 


64 x 64 x 64 

10 

8.14(-9) 

1.00(-4) 

2.00 

6.95(-8) 

4.00 

3.18(—6) 

2.99 

128 x 128 x 128 

18 

7.68(-9) 

2.51 (-5) 

2.00 

4.39(-9) 

3.98 

3.99(-7) 

3.00 

256 x 256 x 256 

3 

6.26(-9) 

6.28(-6) 

2.00 

1.56(-9) 

1.50 

4.99(-8) 

3.00 

512x512x512 

3 

6.00(-9) 

1.57(-6) 

2.00 

1.38(-9) 

0.17 

6.25(-9) 

3.00 


tolerance e while “RRe” denotes the corresponding relative residual. And the same notations 
are used in all the following tables. Table 5.4 and Table 5.6 list the number of iterations, the 
relative residual of the numerical solution on each grid, the Z/> error between the FE solution 
Uh and the exact solution u, the Lo error between the extrapolated solution if 1 , and the exact 
solution u, the Li error between the initial guess W/, and the FE solution Uh, all convergence 
rates, and the ratio r/, defined as 


(5.3) 


II Wh ~ Uh lb 
II Uh - u H2 


which is used to measure how good W /, approximates U/,. While Table 5.5 and Table 5.7 give 
all errors in L norm and the corresponding convergence rates. 


Table 5.6: Errors and convergence rates with e — 10 9 in L 2 norm for Problem 1. 


Mesh 

Iters 

RRe 

Ilf//, - 

wlh 

\\Uh - 

m ||2 

HWft - 

UhWi 

Hi 

Error 

Order 

Error 

Order 

Error 

Order 

32 x 32 x 32 

8 

2.24(-10) 

1.42(-4) 


1.96(-7) 


2.54(-5) 


1.79(-1) 

64 x 64 x 64 

9 

2.68(-10) 

3.55(-5) 

2.00 

1.24(-8) 

3.98 

3.18(—6) 

2.99 

8.96(-2) 

128 x 128 x 128 

16 

7.56(-10) 

8.87(-6) 

2.00 

7.83(-10) 

3.99 

3.99(-7) 

3.00 

4.50(-2) 

256 x 256 x 256 

78 

8.98(-10) 

2.22(-6) 

2.00 

4.98(—11) 

3.97 

4.99(-8) 

3.00 

2.25(-2) 

512x512x512 

3 

7.13(—10) 

5.55(-7) 

2.00 

3.06(—11) 

0.70 

6.25(—9) 

3.00 

1.13(—2) 


As we can see that initial guess Wh is a third order approximation of the FE solution Uh, 
which validates our theoretical analysis in Section 4.4. And the numerical solution Uh reaches 
the full second accuracy for both tolerances, while the extrapolated solution Uh reaches fourth 
order accuracy on the coarse grids and starts to lose accuracy on fine grids. This is due to the 
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Table 5.7: Errors and convergence rates with e — 10 9 in norm for Problem 1. 


Mesh 

Iters 

RRe 

\\Uh - 

w| |oo 

Wh~u 

|oo 

II w h - 

UhWoo 

Error 

Order 

Error 

Order 

Error 

Order 

32 x 32 x 32 

8 

2.24(-10) 

4.02(-4) 


1-11 (—6) 


6.95(-5) 


64 x 64 x 64 

9 

2.68(-10) 

1.00(-4) 

2.00 

6.95(—8) 

4.00 

8.62(-6) 

3.01 

128 x 128 x 128 

16 

7.56(-10) 

2.51(—5) 

2.00 

4.35(-9) 

4.00 

1.07(-6) 

3.01 

256 x 256 x 256 

78 

8.98(-10) 

6.27(-6) 

2.00 

2.88(-10) 

3.92 

1.34(—7) 

3.00 

512x512x512 

3 

7.13(—10) 

1.57(-6) 

2.00 

1.15(—10) 

1.32 

1.67(-8) 

3.00 


fact that the extrapolated solution Ui, is obtained from two second order numerical solutions 
u/, and U 2 h, these two solutions must be extremely accurate in order to obtain a fourth order 
accurate solution Uh- As the grid becomes finer, the tolerance needs to be more cruel. Thus, 
the extrapolated solution Uh starts to lose convergence order when the grid is fine enough 
since a uniform tolerance is used in our ECMG algorithms. In addition, when the uniform 
tolerance becomes smaller, of course, the extrapolated solution Uh reaches the full fourth 
order accuracy on more grids. As shown in Table 5.7, on the grid 256 x 256 x 256 and 
512x512x512 with the tolerance e = 10 -9 , the maximum error between the extrapolated 
solution Uh and the exact solution u reaches 0(10 I0 ). The accuracy of Uh is quite good, 
although the convergence order does not reach the full fourth order on the finest grid. If we 
want to obtain the fourth order accuracy on the finest grid, an even smaller tolerance should 
be used. 

We further point out that the number of iterations is reduced most significantly on the 
finest grid, which is particularly important when solving the large system. Below we provide 
a short illustration. Since we are using the stopping criteria which is related to the relative 
residual error in L 2 norm, in the following, we will show that the ratio /•/, measures how good 
Wh approximates Uh and qualitatively reflects the number of iterations needed in the JCG 
solver with the initial guess Wh in order to get the full second order accurate solution. Since 
||Wh - Uh\\i is third order convergent and \\Uh - w|b is second order convergent, their ratio r/, 
converges linearly to zero as the mesh is refined (see the last column of Table 5.6). Suppose 
the numerical solution on the previous grid (with mesh size 2h) reaches the full second order 
convergence, then we have 

II W h - u || 2 <|| W h - U h lb + II U h - u || 2 , 

< (1 + r h ) || U h - u || 2 , 

(5.4) < II U 2h -u || 2 ■ 

If we use the initial guess Wh directly as the numerical solution on the grid level h without 
any iterations, the convergence order of Wh is given by 

(5.5) order = log 2 >2- log 2 (l + r h ), 

II W h - u lb 

which means that the smaller the r/, is, the more the convergence order 2 - log 2 (l + r/,j 
approaches 2, the less the iterations required to achieve the full second order accuracy are. 
Since r/, decreases half when the grid is refined once, the convergence order 2-log 2 (l+r/ ? ) will 
reach maximum on the finest grid. Therefore, the number of iterations will be reduced most 
significantly on the finest grid. And this is particularly important when solving large linear 












16 


K. J. PAN AND D. D. HE AND H. L. HU 



(a) k = 0 


(b) k = 1 
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o 



Fig. 5.1: U k l/n - U 1/32 on the plane z = 0.5. 


systems. For this example, we have r/, ~ 0.00113 on the finest grid 512x512x512 for both 
tolerances. If using W/, directly as the numerical solution on the finest grid, the convergent 
rate is already greater than (or equal to) 2 - log 2 (l + r/,) = 2-log 2 (l + 0.00113) 1.9838, 

which is almost the full second order convergent. Therefore, only 3 iterations are needed 
to obtain the full second order convergent results, see the last row and second column of 
Table 5.4-5 .7. 

Moreover, Figure 5.1 presents the contour of the error (on the plane z — 0.5) between 
exact FE solution t/ 1/32 and iterative solutions U k ^ 2 ,k — 0,1,3,8 of the JCG solver with 
e = 10” 9 , where C /^ 32 = W 1/32 is the initial guess. Since the error between the VT 1/32 and f/ 1/32 
is symmetric and oscillate, the high frequency error components can be smoothed out easily. 
As we can see that the high-frequency oscillation is smoothed out after only 3 iterations and 
the error has decreased one order of magnitude. In fact, the error between the 8-th iterative 

solution Ufi 32 and C/1/32 is in the order of 0(10 1[ ) while the relative residual ^ 

is 2.24 x 1CT 10 , which is less than 1CT 9 as shown in the first row and third column of Table 5.6. 

Thus, U j/ 32 is actually the numerical solution U 1 / 32 - 
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Table 5.8: Errors and convergence rates with e — 10 12 in L 2 norm for Problem 2. 


Mesh 

Iters 

RRe 

Ilf//, - 

w|b 

Ilf//, - 

u\\i 

II w h - 

Uh\h 

fh 

Error 

Order 

Error 

Order 

Error 

Order 

40 x 16 x 20 

55 

6.84(-13) 

2.97(-4) 


4.81(-6) 


5.93(-4) 


2.00(+0) 

80 x 32 x 40 

81 

8.84(-13) 

7.50(-5) 

1.99 

3.07(-7) 

3.97 

7.44(-5) 

2.99 

9.92(— 1) 

160 x 64 x 80 

137 

9.47(-13) 

1.89(-5) 

1.99 

1.93(-8) 

3.99 

9.33(-6) 

3.00 

4.94(-l) 

320 x 128 x 160 

136 

9.82(—13) 

4.73(-5) 

2.00 

1.97(-9) 

3.30 

1 -17(—6) 

3.00 

2.47(— 1) 

640 x 256 x 320 

12 

9.95(—13) 

1.18(—6) 

2.00 

2.21 (—9) 

-0.17 

1.46(-7) 

3.00 

1.24(—1) 


Table 5.9: Errors and convergence rates with e = 10 12 in L m norm for Problem 2. 


Mesh 

Iters 

RRe 

Ilf//, - 

|cX5 

Ilf//, - 

w||cx) 

l|W A - 

f//,iu 

Error 

Order 

Error 

Order 

Error 

Order 

40 x 16 x 20 

55 

6.84(—13) 

8.06(-4) 


3.50(-5) 


2.23(-3) 


80 x 32 x 40 

81 

8.84(-13) 

2.02(-4) 

2.00 

2.50(-6) 

3.81 

2.78(—4) 

3.01 

160 x 64 x 80 

137 

9.47(-13) 

5.04(-5) 

2.00 

1.68(-7) 

3.90 

3.47(-5) 

3.00 

320 x 128 x 160 

136 

9.82(—13) 

1.26(-5) 

2.00 

1.06(-8) 

3.98 

4.34(-6) 

3.00 

640 x 256 x 320 

12 

9.95(—13) 

3.14C-6) 

2.00 

6.68(-9) 

0.67 

5.43(-7) 

3.00 


Problem 2. The test Problem 2 can be written as 

d 2 u d 2 u d 2 u o 3n n , 

(5.6) di? + Qy 2 + = ^ -2.57r 2 )exp(z)sin(yx)sin(-;y), inQ = [ 0 , 1 ]\ 

where the boundary conditions are 


(5.7) 


j^(i ,y,z) = j-(x, l ,z) = o 

Oil Oil 


and 

(5.8) 

u(0,y,z) = u(x, 0, z) = 0, u(x,y,0) = sin(—x) sin(—y), u(x,y, 1) = e sin( — x) sin( — y). 


The analytic solution of Eq.(5.9) is 


m(x, y, z) = exp(z) sin( y x) sin( y), 

where the exact solution is a sufficiently smooth function which changes more rapidly in the 
x direction than in the y and z directions. 

Since the solution has the fastest change in the x direction and the slowest change in the y 
direction, we use the coarsest grid 10 x 4 x 5 in the ECMG algorithm. Table 5.8 and Table 5.9 
list the numerical data obtained by ECMG/,,, using a tolerance e = 10~ 12 . Again, initial 
guess Wh is a third order approximation of the FE solution f//,. And the numerical solution 
Uu reaches the full second accuracy, while the extrapolated solution Ui, reaches fourth order 
accuracy but starts to loss accuracy on the last two fine grids since we are using a uniform 
tolerance e = 10 12 on each level of grid. Additionally, the ratio r/, converges to zero with 
order one. If using W/, as the numerical solution on the finest grid 640 x 256 x 320, the 
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Table 5.10: Errors and convergence rates with e = 10 1 in L 2 norm for Problem 3. 


Mesh 

Iters 

RRe 

lie,, - 

m|| 2 

lie,, - 

«ll2 

\\W h - 

Uhh 

Th 

Error 

Order 

Error 

Order 

Error 

Order 

32 x 32 x 32 

53 

8.33(—12) 

2.80(-5) 


2.25(-6) 


3.23(—5) 


1.15(+0) 

64 x 64 x 64 

74 

9.95(-12) 

7.16(-6) 

1.97 

2.88(-7) 

2.97 

4.56(-6) 

2.83 

6.37(-l) 

128 x 128 x 128 

52 

9.32(—12) 

1 .8 1(—6) 

1.98 

3.65(-8) 

2.98 

6.25(-7) 

2.87 

3.45(-l) 

256 x 256 x 256 

22 

8.52(—12) 

4.57(-7) 

1.99 

5.14(—9) 

2.83 

8.42(-8) 

2.89 

1.84(—1) 

512x512x512 

9 

7.60(—12) 

1 -16(—7) 

1.98 

2.41(—9) 

1.09 

1.12(—8) 

2.91 

9.66(-2) 


convergent rate is already greater than (or equal to) 2 - log 2 (l + 0.124) ~ 1.8314. Thus, 
only 12 iterations are required to achieve the full second order accuracy, see the last row and 
second column of Table 5.8 and table 5.9. 

Problem 3. Consider a singular solution u e H 3 ~ e (Q) (s is any positive constant) satis¬ 
fying 


(5.9) 


-A u 


33 xyz 

4(x 2 + y 2 + z 2 ) 1 ' 4 ’ 


u = g(x,y,z). 


in O = [0, l] 3 , 
on <90, 


where g(x, y, z) is determined from the exact solution 

u(x,y,z) = — -T-r-TT- 

(x~ + y 2 + z 2 ) 3 ' 4 

The exact solution u has a removing singularity at the origin and has only finite regularity in 

H 2 ~ e . 

Once again, we use 7 level of grids with the coarsest grid 8x8x8 and a tolerance 
e - 10 11 . Table 5.10 lists the numerical data in sense of Li norm starting from the third 
level of grid, i. e., 32 x 32 x 32. The initial guess Wh is still a third order approximation 
of the FE solution (/,„ see second last column in Table 5.10. And the numerical solution (//, 
reaches the full second accuracy. Since the exact solution only has finite regularity in H 2 ~ e , 
the extrapolated solution (7/, can only reach third order accuracy, rather than fourth order 
accuracy for smooth solutions. Moreover, ry, converges linearly to zero. If using Wh as the 
numerical solution on the finest grid 512x512x512, the convergent rate is already greater 
than (or equal to) 2 - log 2 (l + 0.0966) » 1.87. Thus, only 9 iterations are needed to obtain 
the full second order accurate solution, see the last row and second column of Table 5.10. 

Since extrapolation are based on asymptotic error expansions of the FE solution, from 
Table 5.10 it seems that our ECMG method is still effective for such singular problems ( u e 
H ' ~' ; ), and extrapolation can also help us to increase the order of convergence to 3. This is a 
surprising result, which would widen the scope of applicability of our method. 

6. Conclusions. In this paper, we developed a new extrapolation cascadic multigrid 
method, i.e., ECMG 2Cg , for solving the 3D elliptic boundary value problems on rectangular 
domains. The major advantage of our method is to use the Richardson extrapolation and 
tri-quadratic Serendipity interpolation techniques for two numerical solutions on two level of 
grids to obtain a quite good initial guess for the iterative solution on the next finer grid, which 
greatly reduces the iteration numbers for JCG solver. In addition, a relative residual tolerance 
introduced in this paper can be used to control the accuracy of the numerical solutions more 
conveniently, and by using two second order numerical solutions on two scale grids, the 
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fourth order extrapolated solution Uj, on the fine grid can be obtained cheaply and directly. 
Moreover, numerical results show that ECMGj f? has much better efficiency compared to 
classical MG methods and is particularly suitable for solving large scale problems. 

Our method developed in this paper can be easily extended to solve other related equa¬ 
tions, for examples, convection-diffusion equations or Helmholtz equations. Moreover, the 
FE discretization method can be replaced by some other high order methods, such as compact 
finite difference methods. We are currently investigating these extensions. 
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